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Abstract. Pulsars, formed during supernova explosions, are known to be sources of relativistic magnetized winds 
whose interaction with the expanding supernova remnants (SNRs) gives rise to a pulsar wind nebula (PWN). 
We present spherically symmetric relativistic magnetohydrodynamics (RMHD) simulations of the interaction of 
a pulsar wind with the surrounding SNR, both in particle and magnetically dominated regimes. As shown by 
previous simulations, the evolution can be divided in three phases: free expansion, a transient phase characterized 
by the compression and reverberation of the reverse shock, and a final Sedov expansion. The evolution of the 
contact discontinuity between the PWN and the SNR (and consequently of the SNR itself) is almost independent 
of the magnetization of the nebula as long as the total (magnetic plus particle) energy is the same. However, 
a different behaviour of the PWN internal structure is observable during the compression-reverberation phase, 
depending on the degree of magnetization. The simulations were performed using the third order conservative 
scheme by Del Zanna et al. (2003). 

Key words. ISM: supernova remnants - Stars: pulsars: general - Stars: winds, outflows - Magnetohydrodynamics 
- Methods: numerical - Shock Waves - Relativity 



1. Introduction 

Supernova remnants (SNRs) are nebulae originated by the 
explosions of massive stars, known as Supernova events 
(SN), which tipically release an amount of energy ~ 10 53 
erg. Most of this energy is carried away by neutrinos 
produced during the core-collapse phase and the forma- 
tion of a degenerate stellar remnant (neutron star). The 
remaining energy (about 1%) gives rise to a blast wave 
that sweeps up the outer layers of the star and produces 
a strong shock propagating in the surrounding medium. 
The details of such an expansion depend on a number of 
different parameters: the ejected mass and energy, the 
nature and density distribution of the ambient medium 
(jDwarkadas fc Chevalier 19981 IFeatherstone et al. 200T1 
IBlondin et al. 1996|l . the efficiency of radiative losses 
IjUhevalier fc Blondin 1995|) . anisotropies in the supernova 
explosion jWang et. al. 2002| lUhevalier fc Soker 1989|) . 
as well as neutron star spin-down power and proper 
motion space velocity QChatterjee fc Cordes 2002] 
IFrail et al. 19941 IStrom 1 987). Hence, in principle one 
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should expect to find a large variety of different structures 
among PWN-SNR systems. 

If the stellar remnant is a rapidly spinning magnetized 
neutron star (pulsar), then a late energy input is sup- 
plied to the nebular remnant in the form of a wind com- 
posed of relativistic particles (mainly electron-positron 
pairs, with possibly a minority of ions IjAmato et al. 2 003 
and references within)) and a toroidal magnetic field. 
Most of the pulsar rotational energy goes into the launch- 
ing of this believed to be highly relativistic wind, with 
a bulk Lorentz factor tipically estimated to be in the 
range 10 4 - 10 7 iR.ees fc Gunn 19741 IMicnel fc Li 19991) . 
The detailed magnetospheric physics that is at the ori- 
gin of such an outflow is still poorly understood, but one 
point on which all current theories agree is that the wind 
should be magnetically dominated at a distance from the 
pulsar corresponding to the so called light cylinder ra- 
dius, i?LC = c /^j with f2 the pulsar spin frequency. At 
larger distances a bubble of relativistically hot magne- 
tized plasma is then created. This bubble (often called 
"plerion") shines in a very large range of frequencies, from 
radio wavelengths up to X-rays and even 7-rays, due to 
the synchrotron and Inverse Compton emission of the rel- 
ativistic particles. The global properties of the pulsar wind 
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nebula (PWN) can be modeled in terms of injection pa- 
rameters and evolutionary effects l|P acini fc Salvati 19731 
[Reynolds fc Chevalier 1984| ) . However the energy released 
by the pulsar during its entire life-time (~ 10 49 erg in the 
case of the Crab Pulsar) is much less than that driving the 
expansion and evolution of the SNR. Therefore we expect 
the PWN to have limited effects on the global dynamics 
of the SNR. 

The evolution of the PWN-SNR can be divided 
into various phases with different observational features 
flWoltjer 1972| ICioffi 19901 |Reynolds fe Chevalier 1984| 
IChevalier 1998fl . which result from the interaction of the 
SNR structures with the pulsar wind bubble. The SN ex- 
plosion leads to the formation of a SNR in which three 
main discontinuities can be found: a forward shock that 
propagates in the ISM, a reverse shock that propagates in- 
side the free expanding ejecta, and a contact discontinuity 
separating the shocked ISM from the shocked SN material. 
Initially the overall structure expands and also the reverse 
shock moves away from the origin of the explosion, but, 
as the amount of swept ISM grows, the speed of the con- 
tact discontinuity decreases and, finally, the reverse shock 
starts moving back to the origin, where it collapses on a 
time scale of order 10 3 — 10 4 years. Initially (up to about 
10 3 years after the SN explosion, before the reverse shock 
in the ejecta starts moving back), the pulsar spin-down 
luminosity is nearly constant and provides a steady en- 
ergy injection in the PWN. The PWN expands as i 6 / 5 in 
the freely expanding supernova ejecta, in the case where 
the SN has a constant density profile (we will refer to this 
phase as free expansion phase due to the fact that the 
PWN is inside the freely expanding ejecta). The second 
phase (often referred as Sedov- Taylor phase) results when 
the PWN contact discontinuity reaches the SNR reverse 
shock in the ejecta. At this point the PWN expansion 
stops and the SNR reverse shock starts compressing the 
PWN toward the pulsar IMcKee 1974IICiofH et al.T988l) . 
The main observational signature of this process is an in- 
crease in radio luminosity associated to the magnetic field 
enhancement and particle re-energization due to compres- 
sion. In the absence of a PWN the SNR reverse shock 
would collapse to the center, but when a plerion exists 
the presence of an hot bubble prevents the collapse from 
happening and after a transient phase, characterized by 
oscillations of the contact discontinuity, the nebula en- 
ters a phase of slow expansion (Sedov expansion phase 
l|Sedov 1959j> L and finally dissipates in the interstellar 
medium (ISM). 

Hydrodynamical simulations aimed at investigating 
the free expansion and Sedov- Taylor phases were recently 
carried out both in one l|van d er Swal uw et al. 2 001 (SW)) 
and two <|Blondin et al. 20011 (BCF), l.lun 1998)1 dimen- 
sions. These studies were intended to provide some in- 
sights into the details of the system evolution and its sta- 
bility properties (Rayleigh- Taylor instability may occur 
both in the free expansion Ijjun 1998(1 . and in the Sedov- 
Taylor (BCF) stages as a consequence of ejecta acceler- 
ation and compression by the reverse shock). However, 



all the previous investigations dealt with non-relativistic 
hydrodynamic regimes. Our aim is to extend such simu- 
lations to the more realistic relativistic magnetohydrody- 
namical (RMHD hereafter) regime, in order to evaluate 
if the hydrodynamical (HD hereafter) approximation is 
good enough and what differences one may expect. In this 
paper, we present one-dimensional RMHD simulations of 
the first two evolutionary stages in spherical symmetry, 
thus extending the HD study by SW. Two-dimensional 
simulations of the free-expansion phase will be presented 
in a forthcoming paper. We will show that no significant 
differences in the global evolution arise between the HD 
and RMHD cases, both in the free expansion phase and 
during the reverse shock reverberation, even if different 
structures may form inside the PWN. We expect major 
differences to arise in the multidimensional case were the 
toroidal magnetic field may play an essential role in the 
development and growth of instabilities, eventually remov- 
ing degrees of freedom of the system. 

2. Numerical Simulation 

All the simulations have been performed by using 
the newly developed scheme by Del Zanna et al. 
(|I3el Zanna fc Bucciantini 20021 IDel Zanna et al. 2003|l . 
We refer the reader to the cited papers for a detailed 
description of the code, the equations and algorithms 
employed. This is a high resolution conservative (shock- 
capturing) code for 3D-RMHD based on accurate third 
order reconstruction ENO-type algorithms and on an ap- 
proximate Riemann solver flux formula (HLL) which does 
not make use of time-consuming characteristic decompo- 
sition. The code used here has been modified by adding a 
new equation for a tracer that allows us to use different 
adiabatic coefficients for the pulsar wind material and for 
the SNR-ISM medium, namely 7 = 4/3 and 7 = 5/3 re- 
spectively, following the same approach as in Bucciantini 
(2002), but modified for the HLL solver actually used by 
the code. The use of two different adiabatic coefficients 
on a contact discontinuity with a very large density jump 
(density may change by factors over 10 6 — 10 7 , see Fig. 2) 
leads to the formation of waves that tend to propagate 
back into the PWN ( |Shyue 1998] IKarni 1994)l : however 
such waves are well subsonic (the ratio between velocity 
fluctuations and the sound speed is about 0.07 in the HD 
case, and that between kinetic and thermal energy is less 
than 10~ 3 ) so that the behaviour of the PWN is not very 
much affected. 

2.1. Choice of initial parameter and injection 
conditions 

The simulations were performed on a 1024-cell radial grid, 
corresponding to a physical domain extending from the 
origin to 30 Ly. The evolution of the system is followed for 
30,000 years. We set continuous conditions at the outer 
boundary (zeroth order extrapolation) and reflection at 
the origin. Density, momentum-energy and magnetic field 
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from the pulsar are injected in the first cell. No radiation 
cooling is assumed. 

Initial conditions are similar to those adopted by SW: 
total energy in the ejecta E tot — I0 51 erg, mass in the 
ejecta M ej = 3M , ISM mass density lCT 24 g/cm 3 , cho- 
sen to match the partameters for the Crab Nebula. The 
supernova ejecta are set at the initial time in the first 
computational cells: we adopt here a spatially constant 
density (p ) profile and a velocity proportional to the dis- 
tance from the origin such that: 

rRo 

M ej = / 4ir Po r 2 dr (1) 
Jo 

Etot^J o ° ^- Po v 2 r 2 dr (2) 

where R is the initial radius of the ejecta (we have set 
R = 1.2 Ly for a good resolution, but simulations give 
nearly the same results also for smaller values). While in 
the work by SW the supernova energy is initially set as 
enhanced thermal pressure, here we have decided to set 
it as kinetic energy for numerical convenience, in order 
to reduce the initial diffusion of the contact discontinuity 
and thus to obtain a sharper density profile (this problem 
is common to all central-type schemes that avoid spectral 
decomposition) . 

The pulsar wind is created by injecting mass, 
momentum-energy, and a purely toroidal magnetic field in 
the first computational cell, with a total luminosity that 
depends on time in order to include the spin-down process: 

where L D = 5 x I0 38 erg/s, and r = 600 years; no mag- 
netic field is initially present in the SNR or in the ISM. In 
all the simulations we have kept constant in time the fol- 
lowing ratios of injected quantities: the ratio of magnetic 
energy and total energy (a), and that of particle energy 
and mass ~ 100. Units where c = 1 are used. Three dif- 
ferent simulations will be considered here: 

— Purely hydrodynamic; mass and momentum-energy in- 
jected as a wind with Lorentz factor 100, and p/p = 
0.01, no magnetic field. 

— Slightly magnetized wind (a ~ 0.003); as in the pre- 
vious case mass, momentum and particle energy are 
injected as a wind with Lorentz factor 100, and p/p = 
0.01. 

— Magnetically dominated; in this case there is no self 
consistent quasi-steady shock solution on the timescale 
of PWN-SNR evolution (as soon as a shock is formed 
in the PWN, it starts collapsing back to the center). 
We have injected in the first cell magnetic and thermal 
energy with the ratio between the two close to equipar- 
tition (corresponding to a <~ 0.5). In this case we have 
a hot source instead of a cold wind, with p/p ~ 100. 

The injection of the wind is actually a delicate pro- 
cess, especially in the highly magnetized case. The reason 



is that there are cases when the flow in the first cell be- 
comes subsonic, due to the collapse to the center of the 
termination shock (this happens at very early times when 
a highly magnetized wind is considered, but it also hap- 
pens in the low and zero magnetization cases, although at 
later times, since the decreasing pulsar input leads to sit- 
uations in which the shock has to move very close to the 
pulsar in order for the wind ram pressure to contrast the 
inward push from the outer nebula). When the flow in the 
first computational cell is subsonic, there is no longer a 
complete freedom in the choice of the injected quantities. 
We have chosen to assign the values of mass, total energy 
fluxes and a. This forced us to treat the magnetic field not 
as a primitive variable but rather as a derived quantity. 
This is especially important in the magnetic case to en- 
sure energy conservation. Assigning the total energy flux 
and a fixes the magnetic luminosity E mag , so that: 

Emag (t + dt) = E mag (t) + dt X (E mag ) , (4) 

and the value of B is then derived using the following 
equality: 

E mag {t + dt) = 0.5 x {Bit) +dt x Bf. (5) 

We have checked that such injection condition ensures the 
correct energy balance in the PWN. 

3. Discussion 

In this section we will briefly review the various phases 
of the PWN-SNR interaction (the reader is referred to 
SW and BCF for a more detailed analysis in the hydrody- 
namic case), focusing on the effect of the magnetic field, 
especially inside the PWN, which had not been taken into 
account previously. 

First of all, we find an overall slower evolution of the 
SNR, and, consequently, a slower evolution of the PWN, 
with respect to the work by SW, even if all the structures 
appear to be basically the same. As far as the PWN is con- 
cerned, the slower evolution is mainly the effect of a differ- 
ent adiabatic coefficient (4/3 instead of 5/3) that makes 
the nebula more compressible (the PWN size is about 10% 
larger if 5/3 is used, given the same energy input). As for 
the SNR, the different time scale of the evolution is most 
probably due to the different setup of the initial condi- 
tions (the energy released by the explosion is now in the 
form of kinetic energy rather than thermal) . 

We do not expect the global evolution of the system to 
be different between the HD and RMHD case. The only 
two parameters that rule the evolution for a given SNR, 
are the PWN energy and its radius. In SW it is shown that, 
for a given SNR, the evolution of the contact discontinu- 
ity (Rp Wn ) is completely determined by the local value of 
the PWN pressure. When the RMHD case is considered, 
the total pressure may vary inside the PWN due to mag- 
netic tension, but its value at the boundary still stays the 
same since it is a function of E to t alone. Neglecting adia- 
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batic and radiation losses, the following evolution equation 
holds for the PWN in the general RMHD case: 



Density in Log10 — scale, T=1000 years 



Rpwn f{P{Rpwn)) f 



f E tot 



(6) 



where the function / contains the dependences on the SNR 
parameters. 

Let us consider the two extreme cases: no magnetic 
field and magnetic field alone. When no magnetic field is 
present, the total pressure is just the thermal pressure and 
is approximately constant (Ptot = P{Rpwn)) so that: 



UnPtot r 2 dr = 4irP(R pum )R 3 



(7) 



On the other hand, in the magnetically dominated case 
Ptot = S 2 (r)/2, and B(r) ~ 1/r, so: 



4irP tot r 2 dr = 47rP(i? pll)n ) J R^ 



(8) 



and the second equality in Eq. JBJ still holds. 

The thermal and magnetic pressure appear in a com- 
pletely analogous form in the two cases (this is due to the 
fact that the different proportionality coefficient between 
pressure and energy for an HD (1/3) and magnetically 
dominated (1) plasma, exactly compensates the different 
pressure profile) . So it is apparent that what is important 
for determining the propagation speed of the contact dis- 
continuity is just the value of the total pressure, not the 
role that is separately played by the magnetic field and the 
particles. It can be also shown that the relativistic plasma 
and the toroidal magnetic field undergo the same adia- 
batic losses, so the energy variation is analogous for the 
two components under compression or expansion of the 
nebula IjPacini fc Salvati 19 73). As far as radiation losses 
are concerned, these are dynamically important only in 
the pressure dominated portion of the PWN (i.e. only for 
the time a particle spends in the pressure dominated re- 
gion), and not in the magnetically dominated part where 
the dynamics is ruled by the magnetic field. Actually, ra- 
diative cooling may play a more important role in the 
evolution of the SNR. 

Fig. n shows the density and pressure after 1000 
years in the HD simulations: the small oscillations seen 
in the pressure and density profile of the PWN, as 
well as the wall heating (point b) at the boundary 
l|Del Zanna fc Bucciantini 2002|l . are mainly due to the 
use of two different adiabatic coefficients (for a discus- 
sion of spurious wave generation in multifluid simulations 
see IKun fc jishan 1998|) . and they completely disappear 
when a single fluid is used. The various discontinuities 
are all well developed: moving from the origin outward, 
we find the pulsar wind termination shock, the contact 
discontinuity separating the hot relativistic bubble from 
the ejecta material, the shock in the ejecta driven by the 
PWN expansion, the reverse shock propagating into the 
ejecta, the contact discontinuity between the ejected and 
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Fig. 1. Density and pressure profile for the hydrodynamic 
case at time T=l kyr after the supernova explosion. From 
the origin one can see the termination shock of the rela- 
tivistic wind (a), the PWN contact discontinuity (b), the 
swept up shell of ejecta (c), the reverse shock (d), the 
contact discontinuity between ejected and shocked ISM 
material (e) and the blast wave propagating in the unper- 
turbed ISM (f). 



ISM material, and the blast wave propagating in the un- 
perturbed ISM. 

In Fig. the position of the contact discontinuity is 
shown for the HD case. The weakly and strongly mag- 
netized cases are basically the same except for some mi- 
nor differences mainly due to the numerical diffusion of 
the magnetic field. The only way to reduce this effect is 
to use more sophisticated wave-based Riemann solvers. 
However, we do not deem this necessary since the effect 
is kept within a few computational cells. As anticipated 
the behaviour is very similar to that found by SW, even 
if the evolution turns out to be slower. This agrees with 
what one would expect, given the same initial energetics. 
The free expansion phase lasts for about 2000 years, un- 
til the swept up shell of ejecta hits the reverse shock of 
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Evolution of the PWN contact discontinuity 




Density in Log 10 -scale 



Time (years) 

Fig. 2. The radius of the PWN contact discontinuity as 
a function of time: T = 0-1.8 kyr free-expansion phase; 
T= 1.8-15 kyr unstable reverse shock reverberation phase; 
T= 15-30 kyr Sedov expansion phase. The curve shows 
the evolution for the hydrodynamic case, only differences 
too small to be shown on plot are present in the magnetic 
cases, mainly due to magnetic diffusion at the contact dis- 
continuity. 



the SNR. Then the expansion stops and the reverberation 
phase begins: the reverse shock compresses the PWN until 
the pressure of the latter becomes high enough to stop it 
and to push it back. There is a transient phase character- 
ized by oscillations of the PWN contact discontinuity and 
finally the structure relaxes to the Sedov phase. 

In Fig. |3 the time evolution of both density and pres- 
sure is shown. As soon as the reverse shock starts com- 
pressing the PWN, the termination shock moves toward 
the origin and finally collapses to the first computational 
cell. At later times, when the SNR enters the Sedov ex- 
pansion phase, the ram pressure in the wind has dropped 
to such a low value that the termination shock cannot 
be detached from the first cell any more. While the re- 
verse shock moves inward compressing the PWN and in- 
creasing its pressure, the contact discontinuity separating 
the ejecta from the shocked ISM still moves outward. The 
pressure in the ejecta decreases due to rarefaction, and 
when the compression phase has gone on for about half of 
its total duration, it reaches the same values as the pres- 
sure inside the PWN, and the compression starts slowing 
down. As can be seen from the pressure evolution, the 
PWN experienced phases of compression and rarefaction; 
if the latter are strong enough they may give rise to en- 
hancements of the radio emission, so that one can expect 
to see old remnants with high radio luminosity near the 
pulsar position. These oscillations also trigger the forma- 
tion and propagation of weak shocks in the SNR, which 
may reheat the ejecta. These shocks form when compres- 
sions stop and a new expansion phase begins (in our sim- 
ulations we see only two of them, which propagate in the 
SNR with a velocity of the order of 1000 km/s, and a jump 
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Fig. 3. Evolution of a PWN inside a SNR for the HD case. 
Density (top) and pressure (bottom) in logarithmic gray- 
scale and contour levels, with black corresponding to low 
values and white to high values. In the bottom panel the 
position of the contact discontinuity is shown. 



in pressure of about a factor of 2). As shown in BCF the 
number, amplitude and time scale for the oscillations dur- 
ing the reverberation phase, may vary in different PWN- 
SNR systems so that the effect of such weak shocks may 
be important in the case of low luminosity pulsars or light 
SN ejecta. 

In spite of the fact that the evolution of the contact 
discontinuity between the PWN and the ejecta (and as a 
consequence the evolution of the SNR) is the same in both 
the HD and RMHD cases, the structure inside the PWN 
shows different behaviors. This is due to the different stiff- 
ness of a magnetically dominated plasma with respect to 
a purely HD one. For a hot relativistic plasma such differ- 
ences can be easily understood looking at the wave speeds: 
in the magnetically dominated case the wave speed is c 
while in the HD case it cannot exceed c/vo, therefore the 
magnetically dominated region acts more similarly to a 
rigid body in transmitting inward information from the 
PWN boundary. This plays an important role during the 
reverberation phase. In a magnetized case the terminal ve- 
locity inside the PWN tends to a positive asymptotic value 
(that should match the velocity of the contact disconti- 
nuity), which increases as the magnetic pressure becomes 
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Fig. 4. Density and pressure profile for the a = 0.003 
case at time T=l kyr after the supernova explosion. The 
SNR structure is the same as in Fig. 2, while a different 
profile for PWN quantities arises: the termination shock 
is not detached; the total pressure at the PWN contact 
discontinuity is the same as in the HD case. The bottom 
panel shows the total (solid line), the thermal (dotted line) 
and the magnetic (dashed line) pressure. 



more and more important. When the PWN reaches the re- 
verse shock in the SNR, the contact discontinuity cannot 
go any further and starts moving back to the origin. The 
magnetically dominated part (which is in asymptotic con- 
dition) also starts to move back and compress the particle 
dominated region, where mass and energy also increase 
as a consequence of the pulsar injection. A high density, 
and high pressure region is formed near the origin, whose 
dimension is a function of the ratio between the injected 
magnetic and thermal energy during the pulsar life (we 
have assumed it to be constant). In spite of this differ- 
ent behaviour, the magnetic field and particle pressure 
combine in a way to give the same total pressure at the 
boundary in the various cases, as we explained. 
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Fig. 5. Density and pressure profile for the a ~ 0.5 case 
at time T=l kyr after the supernova explosion. In the 
bottom panel different curves represent various pressure 
contributions as in Fig. 4 



Fig. 0] and [3] show density and pressure in the slightly 
magnetized and magnetically dominated case after 1000 
years (the same as in Fig. ^| . As anticipated the position 
of the contact discontinuity and the value of the pres- 
sure at the boundary are the same as the HD case, while 
rather different PWN structures arise. The internal struc- 
ture in Fig. 01 is consistent with the Kennel & Coroniti 
model (Ken nel fc Coroniti 1984|l even if we are not able 
to resolve the termination shock (this should be in the 
first cells but numerical diffusion spreads it to the first 
one). The pressure profile shows a central region domi- 
nated by the thermal pressure and an outer magnetically 
dominated zone: at the contact discontinuity the PWN 
has a ratio between magnetic and thermal pressure ~ 10. 
The wall heating effect is still present even if it is less than 
in the HD case. In Fig. [3] the magnetic pressure is oc r~ 2 
(and almost coincident with the total pressure inside the 
PWN) as expected for a magnetically dominated nebula 
(with a purely toroidal magnetic field), only very near to 
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the origin it reaches values close to equipartition with the 
thermal pressure. 

In Fig.Elthe temporal behaviour of density, total pres- 
sure and thermal pressure is plotted for the case when 
a = 0.003: density and pressure in the SNR are very sim- 
ilar to the hydrodynamic case while inside the PWN they 
show different radial profiles decreasing toward the con- 
tact discontinuity. During the various compression and ex- 
pansion phases due to reverse shock oscillations there are 
enhancements of pressure that can produce a high tem- 
perature and high magnetic field region near the pulsar, 
eventually observable at radio wavelengths. Looking at the 
thermal pressure we can sec the effect of the different rigid- 
ity between the magnetically and thermally dominated 
part of the PWN. As compression starts, the magnetic 
region moves suddenly and starts compressing the ther- 
mally dominated portion near the origin. The material 
starts moving toward the pulsar which still continues to 
inject mass and energy so that these accumulate at small 
radii. Finally, when the system enters the Sedov expansion 
phase, matter and thermal energy start being advected 
away from the pulsar again. This behaviour is completely 
absent in the HD case where density remains constant in- 
side the nebula. This fact is important in that it leads 
to expect different observable signature in old SNR-PWN 
systems depending on the magnetization of the bubble. 
The dimension of the overcompressed region can in prin- 
ciple be used to estimate the magnetic field inside the neb- 
ula and, from this, the magnetization in the pulsar wind 
itself: the larger the size of the overcompressed region the 
lower the magnetic field in the bubble. In the case of no 
magnetic field wc find a uniformly hot bubble, while in a 
magnetized case the emission enhancement may be higher 
near the pulsar and, in principle, even easier to detect. 

A similar behaviour is present in the magnetically 
dominated case shown in Fig.0 but in this case the ther- 
mal pressure becomes too noisy (it must be derived nu- 
merically from the total pressure, which is much larger, so 
that unavoidable accuracy errors arise), and the effect of 
the different rigidities is more easily appreciable by looking 
at the density evolution. The overall nebula can be consid- 
ered to be in a magnetically dominated condition (except 
for the first cell). The injected material remains in the first 
cell during the compression and is advected away only at 
later times when the nebula enters the Sedov expansion 
phase. In this case we expect to find a high pressure and 
high magnetic field spot near the pulsar that may be hot 
enough to be detected not only in radio but also at higher 
frequencies (microwaves or IR). 

3.1. Comparison with existing analytic models 

Two classes of analytic solutions exist in the literature for 
the internal structure of PWNs: the steady-state solution 
IjKennel fc CoronTtr i984 KC) and the self-similar solu- 
tion ( |Emmering fc Chevalier 19871 EC). It is interesting 
to compare the predictions of the analytic models with 
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Fig. 6. Evolution of a PWN inside a SNR for the a = 
0.003 case. Density (top) and total pressure (middle) and 
thermal pressure (bottom) are represented in logarithmic 
gray-scale and contour levels, with black corresponding to 
low values and white to high values. In the middle and 
bottom panels the position of the contact discontinuity is 
shown. 



the results of our simulations, also to evaluate how good 
those models are. We have chosen to compare our results 
with the EC model, since this allows for a non-zero veloc- 
ity of the termination shock. 

It should be emphasized that the applicability of these 
models is limited to the first phase of the evolution of 
the PWN-SNR system, before the reverse shock in the 
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Fig. 7. Evolution of a PWN inside a SNR for the a ~ 0.5 
case. Density (top) and pressure (bottom) are represented 
in logarithmic gray-scale and contour levels, with black 
corresponding to low values and white to high values. In 
the bottom panel the position of the contact discontinuity 
is shown. 



SNR reaches the PWN contact discontinuity. Moreover a 
comparison between our results and the work by EC is 
only possible for the HD and slightly magnetized case, 
the only two case for which a shock solution exists. 

The solution found by EC (as well as that by KC) relies 
on two strong assumptions: a constant pulsar spin-down 
luminosity and a constant velocity at the outer boundary 
of the PWN. None of these applies to our case: in our simu- 
lations the spin-down power fed to the nebula by the PSR 
varies according to Eq.(3) and the velocity of the contact 
discontinuity increases with time during the evolutionary 
stage we are considering, changing from an initial value 
of 0.005c to a value of 0.007c when the system is 2000 
years old. Both these factors lead to differences between 
our results and those found by EC, but surely the different 
energy input as a function of time plays the most impor- 
tant role. As the PSR luminosity decreases, the position of 
the wind termination shock moves to smaller radii than in 
the constant luminosity case. This effect can be easily es- 
timated in the pure HD case. Assuming a constant energy 



input, the position of the termination shock R s (t) can be 
estimated through the condition of pressure equilibrium: 



iwcR 2 s (t) AnR3 pwn ' 

this becomes in the time dependent case: 

Lit) = fiL(t)dt 
A-KcR? s {t) AirR* wn ■ K ' 

When the PSR luminosity is described by Eq.(3), after 
1000 years, the value of R s is a factor 0.612 less than in a 
case when the luminosity is taken to be constant and such 
to give the the same value of R pW n- 

Taking the EC solution in the unmagnetized case and 
with the shock velocity set to 3.5 x 10 _4 c (as evaluated 
from our simulation), we can determine the position of 
the contact discontinuity (normalized to the termination 
shock radius) as the point where the flow speed equals 
the values of 7. x 10~ 3 c, corresponding to the velocity of 
the contact discontinuity in our simulation). Doing this, 
we find (R pwn / 'R s )ec ~ 0-5(Rpwn/Rs)sim, where the sub- 
script ec is the ratio computed on the basis of EC model, 
and the subscript S i m is our simulation value. This is 
consistent with our expectations based on the discussion 
above. 

In Fig.[S]the numerical simulation is compared with the 
EC solution scaled to the same wind termination shock ra- 
dius. Apart from the small waves previously discussed, we 
see a good agreement in the post shock region. The two 
models show bigger differences in the outer part of the 
nebula. As we mentioned, the EC model gives a smaller 
nebular radius: the outer nebula is "more stretched" in 
our case (we want to point out that even if the analytic 
solution exends up to about 5.5 Ly, the dimension of the 
nebula determined by matching the boundary velocity is 
smaller, ~3.6 Ly). This difference is the result of the differ- 
ent pulsar spin-down law: material in the outer part was 
injected at early times, when the luminosity was larger, 
and carries more energy, so that it tends to expand as the 
wind ram pressure drops, and to push the more recently 
injected material close to the pulsar. 

A comparison in the a — 0.003 case is even more del- 
icate: the variation of the boundary velocity is now more 
important because the velocity is close to the asymptotic 
value for a shocked wind with such magnetization so that 
even small variations can lead to major differences in the 
nebular size. Again we have chosen the EC model with 
the same magnetization for a comparison, even if in this 
case the termination shock is not detached from the first 
cell so that we could easily assume that its velocity is 
zero. Results are shown in Figs. 03 and both for the 
overall nebula and for the post shock region. The analytic 
model has been rescaled to the shock position resulting 
from our simulation (determined through a best fit proce- 
dure) which turns out to be at about 0.07 Ly from the ori- 
gin. The radial profiles agree very well with the EC model, 
even if the nebula is more extended and the magnetically 
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Fig. 8. Comparison of the density (upper) and pressure 
(bottom) profiles in the pure hydrodynamic case. The 
solid line is our numerical simulation while the dashed 
line is the EM model rescalcd to the termination shock 
radius. 



dominated part appears to be larger. The differences in 
dimension can again be in part accounted for as an effect 
of the different pulsar spin-down law (i?g should be 0.2 
Ly if L{t) were constant), but likely, in this case, varia- 
tions of the velocity at the boundary as well as numerical 
discretization are also of some importance. 

Our results suggest that, taking into account the spin- 
down process, the Crab Nebula, where the ratio R pwn /R s 
is believed to be about 20, should probably have a smaller 
magnetization parameter than what is estimated by KC. 
It should be noted that EC arrived to the same conclusion 
based on the speed of the termination shock. However this 
subject deserves further investigation. 

4. Conclusion 

In the present paper relativistic MHD simulations of the 
evolution of a PWN inside a SNR, in the spherically sym- 
metric approximation, are shown for the first time. The 
early free expansion and Sedov- Taylor stages have been 
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Fig. 9. Comparison of the density (upper) and thermal 
pressure (middle) and magnetic pressure (bottom) profiles 
for the case a = 0.003. The solid line is our numerical 
simulation while the dashed line is the EM model rescaled 
to the termination shock radius. 
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Same as Fig. 10 but for the immediate post shock 



studied, both in the hydrodynamical and magnetohydro- 
dynamical (with different magnetizations) regimes for the 
PWN. These simulations are mainly aimed at extending 
previous similar works (especially SW) to the more ap- 
propriate relativistic and magnetic regimes. Moreover, two 
different adiabatic coefficients are used for the PWN and 
for the SNR (two fluids). 

An impoertant result is that the evolution of the 
contact discontinuity only depends on the total en- 
ergy in the PWN while it is independent on its mag- 
netization. The same result was previously found by 
|Emmering fc Chevalier 1987| under the assumptions of 
time-constancy of the PSR power input and of the ve- 
locity of the contact discontinuity. 

Our simulations do not take into account radiative 
losses. Such losses may affect the PWN in the early phases 
when synchrotron losses are especially intense, but, in 
principle, on the time scale of the PWN-SNR system evo- 
lution they could be accounted for by renormalizing the 
pulsar spin-down energy input (and eventually changing 
the magnetization of the nebula). Radiative losses may be 
important in a non-trivial way for the dynamics of the 
contact discontinuity and the SNR evolution. Radiation 
cooling may affect the dynamics of the swept up shell of 
ejecta as well as the pressure evolution during the rever- 
beration phase eventually reducing the time-scale and the 
amplitude of the oscillations of the contact discontinuity. 
The main simplification of the model is, however, the re- 
duction of the degrees of freedom of the system because of 
the assumption of spherical symmetry (i.e. 1-D). This pre- 
vents us from addressing the very important issue of how 
the instabilities, that are known to occur at the interface 
between the PWN and the ejecta, both during the expan- 
sion and reverberation phase, develop and evolve. At the 
same time we completely neglect the intrinsic anisotropics 
that the system is likely to posses, like those due to the 
pulsar spatial velocity or to the internal structure of the 
ejecta. 

We also found that the time-dependence of the spin- 
down process of the pulsar and the decrease of ram pres- 
sure in the wind might be an important element in the 
time dependent modelling of PWNs. Steady-state models 
can provide a good description of the radial dependence 
of the various quantities (in units of the distance from the 
shock) but they cannot be used to estimate with confi- 
dence the relative dimension of the contact discontinuity 
and the wind termination shock. This argument however 
deserves further investigation which is beyond the scope 
of this article. 

As in previous studies we find that the transition be- 
tween the free expansion phase and the Sedov- Taylor ex- 
pansion are unsteady and characterized by reverse shock 
reverberation that increases both the pressure and the 
magnetic field in the PWN, eventually reheating the 
plasma and producing bursts of radio emission in the later 
phase of the PWN life, when the pulsar has radiated al- 
most all its spin-down energy. We do not find significant 
differences in the magnetic case. However, magnetic forces 
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inside the nebula give pressure and density profiles with 
higher values near the origin (provided the total energy is 
the same). Moreover, the pressure at the boundary keeps 
the same value and the evolution of the contact discon- 
tinuity and of the SNR is unchanged. Even compression 
does not alter the global behaviour of the system, but the 
presence inside the PWN of two different regions (an in- 
ternal thermal pressure dominated and an outer magnetic 
dominated zone) with different associated rigidities gives 
rise to an overcompression near the origin, where mass 
and energy coming from the pulsar are confined inside a 
small radius. In this case a high-temperature and high- 
magnetic field spot is formed, which could eventually be 
revealed not only in radio but also at higher frequencies 
(microwaves or IR). 
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